//
// Esempio con matrice di Vandermonde. La matrice ha un numero di
// condizionamento rapidamente crescente con la dimensione del sistema.
//
// Il sistema è risolto in due diversi modi: utilizzando EGPP ed utilizzando
// qr.
//
// Osservazione: la misura relativa dei residui N_RES (paragonabile nei due casi)
// indica che il sistema è stato risolto "bene" dal calcolatore. L'errore
// relativo, rapidamente crescente con la dimensione, è dovuto al pessimo
// condizionamento del sistema. 
//
clear;
//
Percorso = "PERCORSO files seguenti";  // <------ *** MODIFICARE ***
exec(Percorso + "EGPP.sci");
exec(Percorso + "SA.sci");
exec(Percorso + "SI.sci");
//
// variabili utili
//
DIM = [];
COND_A = [];
ERR_EG = [];
N_RES_EG = [];
ERR_QR = [];
N_RES_QR = [];
SUP_EG = [];
SUP_QR = [];
//
// ciclo sulla dimensione del sistema
//
for N = 5:2:25,
DIM($+1) = N;
//
// generazione della matrice
//
colonna = linspace(0,1,N)';
A = zeros(N,N);
for k=1:N,
    A(:,k) = colonna.^(k-1);
end;
//
// soluzione esatta e termine noto corrispondente
//
x_esatta = ones(N,1)/sqrt(2);
b = A*x_esatta;
//
// numero di condizionamento della matrice
//
COND_A($+1) = cond(A);
//
// soluzione del sistema con EGPP, SA e SI
//
[S,D,P,info_EGPP] = EGPP(A);
[c_EG,info_SA] = SA(S,P*b,"disabilita");
[x_EG,info_SI] = SI(D,c_EG,"disabilita");
//
// soluzione con qr, moltiplicazione e SI
//
[U,T] = qr(A);
c_QR = U' * b;
[x_QR,info_SI] = SI(T,c_QR,"disabilita");
//
// misura relativa degli errori
//
ERR_EG($+1) = norm(x_EG - x_esatta)/norm(x_esatta);
ERR_QR($+1) = norm(x_QR - x_esatta)/norm(x_esatta);
//
// misura relativa dei residui
//
N_RES_EG($+1) = norm(A*x_EG - b)/norm(b);
N_RES_QR($+1) = norm(A*x_QR - b)/norm(b);
//
// uso del Teorema di condizionamento
//
SUP_EG($+1) = COND_A($)*N_RES_EG($);
SUP_QR($+1) = COND_A($)*N_RES_QR($);
end;
//
// grafici
//
scf(0);clf(0);
plot(DIM,log10(COND_A),"gs",...
     DIM,log10(SUP_EG),"rs",DIM,log10(ERR_EG),"r+",DIM,log10(N_RES_EG),"rx",...
     DIM,log10(SUP_QR),"bo",DIM,log10(ERR_QR),"b+",DIM,log10(N_RES_QR),"bx");
xlabel("dimensione");
ylabel("logaritmo in base 10 di...");
legend("COND_A","SUP_EG","ERR_EG","N_RES_EG","SUP_QR","ERR_QR","N_RES_QR",2);
title("confronto EGPP , qr");
